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SUMMARY 

This report describes the construction of an explicit, single time-step, conservative, finite- 
volume method for multidimensional advective flow, based on a uniformly third-order poly- 
nomial interpolation algorithm (UTOPIA). Particular attention is paid to the problem of 
flow-to-grid angle-dependent, anisotropic distortion typical of one-dimensional schemes used 
component-wise. The third-order multidimensional scheme automatically includes certain 
cross-difference terms that guarantee good isotropy (and stability). However, above first- 
order, polynomial-based advection schemes do not preserve positivity (the multidimensional 
analogue of monotonicity). For this reason, a multidimensional generalisation of the first 
author’s universal flux-limiter is sought. This is a very challenging problem. A simple flux- 
limiter can be found; but this introduces strong anisotropic distortion. A more sophisticated 
technique, limiting part of the flux and then restoring the isotropy-maintaining cross-terms 
afterwards, gives more satisfactory results. Test cases are confined to two dimensions; 
three-dimensional extensions are briefly discussed. 


INTRODUCTION 

Many numerical models in computational fluid dynamics employ so-called “leap-frog” 
finite-difference schemes using second-order centred differences in both time and space. 
One significant attribute of such methods is that for purely advective flow there is no 
artificial dissipation supplied from the numerical scheme. This has both positive and 
negative consequences. On the one hand, where the flow is well resolved, the numerical 
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method does not corrupt modelled physical diffusivity or viscosity. However, it is well 
known that in the presence of sharp changes in gradient, centred schemes exhibit serious 
phase errors generating spurious oscillations and unphysical maxima and minima in the 
transported variables. And because of the lack of inherent numerical dissipation, these 
oscillations are not strongly damped and may lead to unnatural processes such as “numerical 
turbulence” thereby producing anomalous mixing. First-order methods are nonoscillatory, 
but contain too much inherent artificial diffusivity to be seriously considered as the basis for 
a multidimensional advection scheme. By contrast, third-order upwind schemes possess the 
inherent stability of upwinding in addition to high accuracy and excellent phase behaviour, 
although small unphysical undershoots or overshoots may still be excited in sharp-gradient 
regions. The current work focuses on the development of an explicit, single-time-step, 
conservative, flux-based, control-volume algorithm for multidimensional advection. This 
is in contrast to the common practice of operator-splitting (using locally one-dimensional 
algorithms) which, by its nature, cannot be formulated in a conservative manner. 

One other aspect of multidimensional advection schemes that has received little 
attention is the problem of anisotropic distortion that occurs when advection is oblique or 
skew to the grid lines. The problem stems from using one-dimensional advection formulae 
in each coordinate direction separately, without explicitly taking account of necessary cross- 
difference terms. In practical simulations of physical processes, the effects of this distortion 
are not generally immediately apparent because the flow patterns are very complex and the 
true solution is not normally known. However, simple model test problems (with known 
solutions) of advection oblique to the grid suggest that serious quantitative errors may be 
introduced through this mechanism in some circumstances. It should be noted that the 
commonly used test problem of advection of an initial profile around a circular path does not 
highlight this type of distortion because the flow-to-grid angle is continuously changing. But 
in a relatively constant-direction flow field, this problem can evidently lead to quite 
substantial errors. In addition, above first-order, one-dimensional schemes used component- 
wise for multidimensional pure advection are unconditionally unstable (unless the flow is 
aligned with the grid). This is usually a weak long-wavelength instability and is often not 
excited in spatially varying velocity fields. Although a small amount of physical diffusion 
is usually enough to suppress the instability, an advectively stable algorithm seems more 
desirable. 

A number of one-dimensional monotonicity-preserving (MP) advection schemes have 
recently been compared (Bull, 1990). Schemes tested include a second-order total-variation- 
diminishing (TVD) flux-limited scheme due to van Leer (1974), referred to as the curved- 
line advection method (CLAM), and Leonard’s (1991) ULTIMATE QUICKEST algorithm, 
a third-order one-dimensional MP advection scheme. These methods have also been tested 
in two dimensions with explicit forward time-stepping using fluxes based on the one- 
dimensional formulae in each coordinate direction. Although the CLAM scheme maintains 
“positivity” (the multidimensional analogue of monotonicity) provided the Courant number 
is small enough, the one-dimensional ULTIMATE scheme, applied coordinate-wise does not. 
In addition, both methods generate strong anisotropic distortion, although of a different type 
to that observed with unlimited centred differences. 

The following sections address the two outstanding problems of (i): eliminating 
anisotropic distortion and (ii): developing a truly multidimensional flux-limiting strategy. 
The proposed techniques are tested on a conceptually simple but numerically very 
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challenging problem involving pure two-dimensional advection of an initial circularly 
symmetric Gaussian profile by a constant velocity field oblique to a uniform square mesh. 
An angle of 45° usually gives worst-case distortion; however, other angles are also 
considered in order to develop some perspective on flow-to-grid angle dependence. All of 
the test cases and most of the discussion will be confined to two dimensions. Although there 
is a large theoretical gap to be bridged between one and two dimensions, it appears that, 
once this is achieved, algorithmic extension to three dimensional advection may be relatively 
s traig h t- forward . 


ELIMINATION OF ANISOTROPIC DISTORTION 

Anisotropic distortion appears in multidimensional advection schemes that rely on the 
respective one-dimensional formulae in each coordinate direction separately. As a result of 
this, computed solutions are strongly dependent on the flow-to-grid angle; equivalently, this 
means that solutions are not grid-orientation independent. This problem can be seen 
theoretically by making use of the two-dimensional complex amplitude ratio (or amplification 
factor) G. For reference, consider a model two-dimensional advection-diffusion equation 
for a scalar <j>(x,y,t) 


= p(d 2 4> ^ d 2 <t>\ 
dt dx dy \ 3jc 2 dy 2 ' 


( 1 ) 


where u, v, and D are constant coefficients. It is of interest to follow the evolution of a 
wave of the form 


4> = A (/) exp(i kx) exp(t k\) 

x y + 


(2) 


involving component wave-numbers k x and k y , where A (t) is the (complex) amplitude and 
t (iota) represents the imaginary unit, V-l . The complex amplitude ratio G is then defined 
as 


G = </>(*, y > f + Ar) 

y, t) 


A(r+Ar) 

W) 


(3) 


Substitution of (2) into (1) gives 

4> = A(0) exp[ - (k x + ky)Dt] exp[i k x (x - ut )] exp[t k y (y - vt )] ( 4 ) 

so that the exact expression for the two-dimensional G of the model advection-diffusion 
equation is 

G «.ct = ex p[-(a x d 2 x ->-a y e 2 y )] exp(-ic^) exp(-i c y 6 y ) ( 5 ) 

where the diffusion parameters 


a 


X 


DAJ 

Ax 2 


( 6 ) 
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and 


and component Courant numbers 


and 


a.. 


D Ar 
Ay 7 


uAt 

Ax 


vAt 

~Ay 


have been introduced, and the nondimensional wave-numbers are given by 


6 X = k x Ax 
Q y = k y Ay 


(7) 

( 8 ) 

(9) 

( 10 ) 

( 11 ) 


For purely advective flow in two dimensions, D = 0 and (5) becomes 

G .dv = ex P["‘(cA + c/y)] (12) 


A Taylor expansion of real and imaginary parts of this expression gives, to third order in 6 X 
and 6 y , 


adv 


= 1 - 


i^e 2 x + c x e x c y e y + ic y d 2 y 


- I 


( C A + C y 0 y) - i A 0 , + 


(13) 


The corresponding numerical G for a given computational method (involving a single 
forward time-step) can be similarly expanded in a Taylor series. It is not difficult to see that 
the use of only one-dimensional advective formulae in each coordinate direction would result 
in a G-expansion involving powers of ( c x 6 x ) and ( c y 6 y ) separately, but no cross-terms. In 
particular, the important second-order term 

M,)«W < 14 > 

appearing in (13), is absent from such schemes. Thus, flow along coordinate directions 
(either c x = 0 or c y =0) may be well modelled, but at other angles the missing term will 
clearly lead to errors. One might expect most distortion to occur for flow along grid- 
diagonals (i.e., at 45° on a uniform square mesh), and this is in fact borne out by numerical 
experiments. 

As will be shown, it is possible to devise an explicit, forward-time-step, conservative 
control-volume formulation of advective transport that includes all the terms (through third 
order) involved in (13). It is also a straight-forward matter to generalise the method to three 
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dimensions (see Appendix 1), and to extend the theory to include diffusion terms to the same 
order (Appendix 2). This uniformly third-order multidimensional scheme is not, however, 
positivity preserving. Rather, it is the analogue of the third-order one-dimensional 
QUICKEST scheme (Leonard, 1979), which also generates (albeit relatively small) over- 
shoots or undershoots when modelling sharp changes in gradient. In the one-dimensional 
case, a rather simple universal limiter (valid for any order of accuracy) can be applied, 
resulting in the monotonicity-preserving ULTIMATE QUICKEST scheme (Leonard, 1991). 
The extension of the limiting strategy to two (and three) dimensions is by no means straight- 
forward; however, some progress has been made, and this is described in the next section. 
The remainder of this section is devoted to studying anisotropic distortion of unlimited 
schemes and to the description of the two-dimensional uniformly third-order advection 
method. 

A simple test problem has been set up to study anisotropic distortion under purely 
advective conditions. This consists of a uniform square grid (Ax = A y = h) over a 31x31 
domain (i.e., h = 1/31) with a uniform velocity field at an angle 0 to the x-axis 

V = ( u , v) = (V o cos0, V o sin0) (15) 

The problem simulates the equation 


^ + y -v<t> = o 

dt 


(16) 


representing pure advection of a simple scalar “density” profile, taken here to be 

a circularly symmetric Gaussian distribution with a standard deviation of 3 h, initially 
established at the centre of the grid. The domain has doubly periodic boundary conditions 
so that the profile can be advected any distance, depending on the run time. In most cases, 
0 is chosen to be 45° ( c x = O, and At is adjusted so that the (exact-solution) profile would 
be advected a distance 31V2 h; i.e., to the centre of the next diagonal “box”. This is 
shown schematically for the exact solution in Figure 1(a), giving contours of constant 
density. Figure 1(b) shows a “three-dimensional” portrayal of the initial (s= final) distri- 
bution. Of course, the actual computation takes place within a single box, with the profile 
“exiting” at the top-right corner and “re-entering” at the lower-left, with the final position 
at the center again. For the 45° case, c x = c y = 31 IN, where N is the number of time steps. 

By contrast to the exact solution of Figure 1, Figure 2 shows the results of using the 
(“standard”) space-time-centred leap-frog scheme, 


i n + 1 

0i, j = 




- ~ +UJ ~ Cyttb* 1 ' 4>lj- 1 ) 


(17) 


for c x = c y = 0.25 (N = 124). Clearly, even for this relatively smooth initial profile, there 
has been significant anisotropic distortion in addition to serious (phase-lag) oscillations due 
to the inherent (undamped) numerical dispersion, characteristic of this scheme. The second- 
order (three time-level) leap-frog scheme gives results showing similar distortion to that of 
a method using the (two time-level) explicit one-dimensional Lax-Wendroff (1960) scheme 
in each coordinate direction separately, 
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4>u = Ki ~ ~ ti-hj) - j(4>U,j - 2*1 j + 4>U,j) 


- C 


- k,-<) - - w.i * *;,/-■) 


( 18 ) 


This scheme is unconditionally unstable unless the flow is aligned with one of the coordinate 
directions. However, for small time-steps, the instability is weak, and results can be 
obtained over short run-times, as shown in Figure 3. Again, significant distortion and phase- 
lag oscillations are immediately obvious. For extended run-times, the computation will 
ultimately “blow-up”. 

As in the one-dimensional case, oscillations can be avoided by reverting to a first-order 
upwind technique. If one employs a simple-minded, coordinate-wise one-dimensional 
scheme (for c x and c y both positive) 


K} = 4>lj ~ c x (<f>l j - 4>U,j) ~ cA*lj - 


(19) 


the results are as shown in Figure 4 for c x = c y = 0.5. In this case there are no undershoots 
(this is the simplest two-dimensional positivity-preserving advection scheme). However, 
gross anisotropic distortion is observed, as is a massive loss in peak value. Both of these 
effects are due to the large amount of artificial numerical diffusion inherent in this method. 
Figure 5 shows the diamond-shaped von Neumann stability region in the ( c x , c y ) plane for 
first-order upwinding used coordinate-wise, allowing for positive or negative velocity 
components; i.e., |c,| + |c y | < 1. 

Equation (19) can be rewritten as a conservative, control-volume, explicit update 
equation 

C 1 = 4>lj + FLUXW(/, j) - FLUXW(i+l, j) 

+ FLUXS(i, j) - FLUXS(/, ; + l) ( 20 ) 


by identifying the west-face advective flux (for c x > 0) as 

FLUXW (i,j) = c x 4> w = c x *U,j 
and the south-face flux (for c y > 0) as 

FLUXS (i,j) = c y 4> s = c y *lj. x 


( 21 ) 

( 22 ) 


where * w and <t> s are considered to be the effective west and south face-values, respectively, 
advected into a given control volume. Note that (20) guarantees conservation by expressing 
the east-face flux in terms of the west-face flux at control-volume (i+l,j), and similarly for 
the north-face flux in terms of the south-face of CV(/, j+ 1). This simple type of first-order 
upwinding assumes that the scalar value advected across any particular face is effectively 
equal to the immediate upwind node-value in a direction normal to that face. But this does 
not take account of the fact that, in general, the advective flux is sweeping across the face 
at an angle (rather than simply in the normal direction). Figure 6 suggests that the (upwind- 
biased) transverse gradient of * should be important in this respect. In fact, if one assumes 
that in the vicinity of the west face, the scalar value varies linearly between and 
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it is not difficult to show that the corresponding average face-value (over a time-step 

At) is 


<t> 


w 




- <t>U,j-i) 


(23) 


thus differing from the simple-minded upwind formula by the appearance of the transverse- 
gradient term (for c x , c y > 0) 

GRADTW = 4>U,j - (24) 

Similarly, the south-face value can be shown to be 

4> s = 4>" j_ x - GRADTS (25) 

where (again for c x , c y > 0) 

GRADTS = (26) 

When these formulae are used in (20), the results are as shown in Figure 7. By comparison 
with Figure 4, it is immediately obvious that inclusion of the transverse-gradient terms has 
restored isotropy (at least for 0 = 45°), although there is still significant artificial diffusion 
as evidenced by the strong loss in peak value and concomitant spreading. At other advection 
angles there is some anisotropic distortion, but not nearly as much as with the simple first- 
order scheme of (19). The von Neumann stability region in the ( c x , c y ) plane is now a 
square: (l-|c Jt |)(l-|c y |) > 0. Also note that exact point-to-point transfer occurs (across 
a grid diagonal) when \c x \ = | c y | = 1 . 

Inclusion of the GRADT terms can be directly related to the appearance of the second- 
order cross-term (c x 6 x )(c y 6 y ) in the Taylor expansion of the numerical complex amplitude 
ratio. Thus it would appear that the two-dimensional, conservative, flux-based, control- 
volume formulation of the Lax-Wendroff scheme would benefit from inclusion of these 
terms, as well. This is suggested in Figure 8, where it is assumed that <t> has a linear 
dependence on both x and y in the vicinity of the west face, collocated (at time-level n) at 
and (for c x and c y both positive). Clearly from the figure, these are the 

appropriate (upwind-biased) nodes. Taking the time average (over At) of the west-face value 
(equivalent to averaging 4> n over the area covered by the small angled arrows) then leads to 


j (K, + j - tU.i) - <27) 


where one again sees the appearance of the GRADTW term, defined in (24). The corre- 
sponding south-face value (for c x and c y both positive) is 


4> s - ^ (4>l j + 4>1,j- 1) - y (<t>l j - 4>1j- 1 ) - y j-i - 4>1-i,j-i) 


(28) 


It should be clear that the GRADT terms are dependent on the signs of both c x and c y , 
although the other (face-centred) terms in (27) and (28) are independent of the signs of the 
respective advecting velocity components. The results of using (27) and (28) in (20) are 
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shown in Figure 9 for c x = c y = 0.25. This is clearly an improvement over Figure 3, 
although there is still a large degree of phase-lag dispersion (just as in the one-dimensional 
case). The Taylor expansion of the numerical G corresponding to this scheme matches (13) 
identically through all second- order terms (including the c x d x c y 6 y cross-term). Thus it can 
be considered to be the two-dimensional, conservative, control-volume analogue of the Lax- 
Wendroff scheme. Inclusion of the transverse-gradient terms also restores stability; the 
stable region is once again the square: (1 - 1 c x | ) (1 - 1 c y | ) >0. Extension to three dimensions 
is straight-forward, simply involving an additional transverse-gradient term (along the extra 
coordinate direction) in the formulae for each of the three face- values. Incorporation of 
standard second-order differences for diffusive fluxes is consistent with the Taylor expansion 
of (5) involving all second-order combinations of the a’s and c’s. 

Flux-Based Formulation 

It is instructive to re-derive the uniformly second-order advection scheme from a rather 
different view-point. This will suggest a way of generalising to uniformly third (and 
potentially higher) order methods. The idea is based on the fact that for pure advection at 
a constant velocity V, the scalar value at a point x at the new time-level is equal to that at 
position x-VAr evaluated at a time At earlier: 

t+At) = 4>(x-\At, t) (29) 

This means that <f> can be predicted at any given (nodal) point by projecting back along an 
advective characteristic and interpolating spatially to evaluate 4> at the origination point at 
the earlier time-level. The accuracy of the scheme then depends entirely on the form of the 
spatial interpolation used for evaluation of the right-hand side of (29). Although the 
algorithm is derived on the basis of a constant V, the update algorithm is subsequently 
rewritten in a flux-based, conservative control- volume form, wherein the individual 
advecting face velocities are allowed to be spatially varying in conformity with a solenoidal 
velocity field. 

Referring to Figure 10, assume that at the earlier time-level 

*(*, y, i ) - *’«, v) (30) 

where £ and rj are local x- and y-coordinates centered in the control volume. In the case 
shown (c x , c y > 0), assume that 

*’((,<!) - C, ♦ C 2 f ♦ C, f 2 ♦ C, V + C s - C s ( V (31) 

and evaluate the C’s by collocating at the six node- values. Note that the inclusion of the 
(i-l, j- l)-node gives the appropriate upwind bias in the quadrant corresponding to the 
position of the origination point (shown by the circled star). It Should be clear how the 
stencil changes for origination points in other quadrants. 
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It is not difficult to show that the C’s are given by 

C t = 4>lj 


(32) 


C 2 = 


(4>Iuj ~ 


2 Ax 


(33) 


^3 = 


(4>luj ~ 2 Kj + 4>hJ 


C 4 = 


2 Ax' 1 

i - Kj- i) 


^5 = 


2 Ay 

~ 20" j + 0" 7 -,) 


and 


Q = 


2 Ay 2 

(0"j - *U,j - tij-i + _ 
Ax Av 


(34) 


(35) 


(36) 


(37) 


The upwind bias of the latter coefficient, of course, depends on the signs of c x and c y , 
whereas the other coefficients are centred with respect to node ( i,j ). It is easy to see how 
C 6 would change for the three other combinations of signs of c x and c y . 

From (29) and (30), the explicit update formula can be written in general as 


= 4> n (-uAt, -vA t) 


(38) 


For the particular (second-order) case shown in Figure 10, substituting the above C-values 
results in 


, n + 1 



*h.j) + y(tf*i./ 

c 

C 2 

^(0",-> - < 

*>",-l) + -f (*?./.! 

Wy&lj ~ * 

, n , n , j 


* n i n 


>-i) (39) 

This can be rewritten in conservative control-volume form, as given by (20), by identifying 
west and south face-values as, respectively, 

- !(*?,; + *?-,./) - j(*1.j - *1-1. i) - ^GRADTW 


(40) 
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and 


*, - * Kj-d - - *1.1- 1) - yGRADTS (“1) 

where GRADTW and GRADTS are given (for c x , c y > 0) by (24) and (26), respectively. 
These face-values are identical to those given by (27) and (28)— i.e., the one-dimensional 
Lax-Wendroff formulae together with the appropriate (upwind-biased) GRADT terms. Note 
particularly how the cross-difference term, given by the last term in (39), has been split 
between x- and y-fluxes. This is made possible because of the identity 


(*<.7 ~ *i-i, j) " (*<, j-i ~ 


(*(, j ~ $i,}- 1) " j ~ 


corresponding to the cross-derivative identity 


m cfy 

By Bx Bx By 


(42) 


(43) 


This type of manipulation is particularly important in developing higher-order formulae, such 
as the uniformly third-order method considered next. It should also be noted from (39) that 
exact nodal point-to-point transfer occurs when c x = 1 and c y = 0 (or vice versa), and also 
(across the grid diagonal) when c x = c y = 1, when GRADT-terms are included. 


Uniformly Third Order Polynomial Interpolation Algorithm 

Figure 11 shows the appropriate upwind-biased 10-point stencil for collocating a third- 
order polynomial of the form 

= c, + c 2 * + c 3 e + c 4 e 

+ C 5 V + c 6 I rj + C 7 I 2 V 

* C 8 V 2 * C 9 Ir I 2 

+ C l0 V 3 (44) 

For notational convenience the subscripts on the node-values denote central -point (P), west 
(W), east (E), north (N), south (S), north-west (NW), double-west (WW), etc. (as usual, 
capital letters are used for node- values and lower-case for face- values); the time-level n 
superscript has also been dropped. In Figure 11, the collocation points have been chosen 
so as to give some upwind bias corresponding to the location of the origination point. Also 
note the symmetry of the stencil about a diagonal closest in direction to that of the advecting 
velocity vector. (It should be clear how the stencil changes for location of the origination 
point in other quadrants). For this choice of stencil, it is not difficult to derive the following 
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C- values: 



c, = <t> p 



(45) 

'O- 

1 

Ul 

ll 

( <t> E ~ 3 <f> p 

+ 3 4> w 

- <t>ww) 

(46) 

2 Ax 


6 Ax 


c, = 

( <t> E - 2<t>p + 

K) 


(47) 

3 

2 Ax 2 




<T> 

II 

X-N 

•e- 

C*1 

- 2>4> P + 3<£ w 
6Ax 3 

- <t>ww) 

(48) 

(<t>N ~ ^s) 

(4> n ~ 3</>, 

. * Its 

- * ss ) 

(49) 

2 Ay 


6 A y 



and 


C 6 = 


(<^£ 4>p ) (&SE + (^tf ^Aw) (^V ^v) 

2 Ax Ay 


(50) 


C 7 = 


(^£ ^Sf) 2(</>p </> s ) + (</> w ^sw) 

2 Ax 2 Ay 


(51) 


C 8 = 


(<£* _ 2 </> p + <t> S ) 

2 Ay 2 


(52) 


C 9 = 


Wv - <A W ) “ 2 (^p " ^w) + (^s ~ ^5iv) 

2 Ax Ay 2 


(53) 


'10 


(<£„ - y>p + 30 5 - <fr 55 ) 

6 Ay 3 


(54) 
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Substitution of these values into (44) and (38) results in the following formula for the 
updated (superscript +) central -point node- value: 


C C x (1 - cl) 


K - K- - K) + 


(<t> E ~ 3</>„ + 3<f> w - 0™,) 


-■f (**-*,) + 7 


C v (l " Cy ) 

l (<t>N ~ 3< £ P + 3< ^5 “ ^Ss) 


+ -y (4> e - 2<t>p + 4> w ) + -y (0* - 2< Ap + <As) 


c c 


+ “V 2 " ^SSE " *S> + ~ ^W) ■ - 4> w )] 


c c 2 


[(^AT " <t>Nw) ~ 2 (<t> P ~ <t> W ) + (4>S “ ^5W>] 


c 

X y 


[(+, - +*) - 2 (*, - *,) * «W - *„)] 


(55) 


This can be written in conservative control-volume form by identifying the effective advected 
west-face value as 


, 1 it X \ ^X i X X \ ^ ^X ) 

tw = 2 ^ + “ y “ 5 


(<t> p - 2 <j> w + 0^) 


[(<£/> <^> s ) + WVw ^ (V )] + y (*^AW 2< ^W + ^sw) 


“T^ ( <t> P ~ <t> W ~ <t> s + <I>SW ) 


(56) 


with an analogous formula for <f> s . Note that the first three terms on the right of (56) 
represent the one-dimensional contribution, corresponding to the QUICKEST scheme 
(Leonard, 1979). The next term can be rewritten using the following identity 


(0p $$) + ($AW (&W 4* SW ) (^AW 2< ^W + ^ SW^ 

4 _ 2 4 

(4> P - <i> w - 4> s + <t> sw ) 


4 


(57) 
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in which case, the latter two terms on the right can be combined with the corresponding 
terms in (56), giving 

c (1 - cl) 

<j> w = 4 > lin - — GRADNW - - — CURVNW 

2 6 

- 5; GRADTW 
2 

- ^[(1 - c y ) CURVTW + (1 - c x ) TWISTW] 5g) 

where, as seen from (56), <t> LIN represents linear interpolation across the west face, 
GRADNW is the normal- gradient term, CURVNW represents the upwind-biased normal- 
curvature term (depending on the sign of c x ), GRADTW is the transverse-gradient term 
appearing earlier in lower-order formulae (depending on the signs of both c x and c y ), 
CURVTW is the upwind-biased transverse-curvature term (depending on the sign of c x ), and 
TWISTW is a cross-difference or “twist” term (again depending on the signs of both 
velocity components at the face). Of course, the south-face value has an entirely analogous 
formula that can be easily derived by simple permutation of node-values. 

Figure 12 shows the results of the standard test-problem using this uniformly third- 
order polynomial interpolation algorithm (UTOPIA). Clearly, there has been a significant 
improvement in performance, as compared with lower-order schemes. Isotropy is very good 
and peak resolution is reasonable, the only significant problem being the appearance of 
(relatively small) negative values in the out-skirts of the main profile region. In some cases, 
undershoots of this kind are merely an inconvenience. However, some transported variables 
may need to be kept strictly non-negative for both physical and (more importantly) stability 
reasons. In such cases, lack of positivity preservation could cause serious problems. 
UTOPIA’S behaviour is relatively insensitive to stream-to-grid angle, as seen in Figure 13 
(0 = tan' 1 0.5 ~ 26.565°) and Figure 14 (0 = 0°). Perhaps surprisingly, the stability 
region is restricted to the diamond-shaped figure: |cj + |c y | < 1 , even though exact point- 
to-point transfer again occurs at |cj = J | = 1 (evidently, these are isolated stable points). 

For comparison, Figure 15 shows 45°-advection results using the (theoretically 
unstable) one-dimensional QUICKEST scheme in each coordinate direction separately; i.e., 
dropping the GRADT-, CURVT-, and TWIST-terms in (58). Such a two-dimensional advec- 
tion scheme, together with standard second-order central differencing for diffusion, has been 
used for many years by Davis and Moore (1982, 1985) and co-workers. Once again, one 
sees considerable anisotropic distortion, although phase error is lower than that of second- 
order methods. For pure advection, this method is also theoretically unstable, unless the 
flow is exactly aligned with the grid. It is perhaps of interest to look at a “quasi-third- 
order” scheme in which face-values are given by the one-dimensional QUICKEST formulae, 
together with only the GRADT terms (i.e., omitting transverse curvature and twist terms, 
in order to effect some economisation). The 45°-advection results for this scheme are shown 
in Figure 16. Although much better than the uniformly second-order method, Figure 9 (in 
terms of phase error), or Davis and Moore’s method (in terms of grid-orientation distortion), 
it is not immediately clear whether or not the savings incurred by algorithmic simplification 
might be more than offset by deterioration of the simulation. In other words, to achieve a 
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given accuracy it might be necessary to refine the grid using the “quasi-third-order” scheme 
to a point where it is actually more expensive, overall, than using UTOPIA on the original 
grid. This question can be answered by performing a series of grid-refinement studies and 
comparing the relative computational efficiency of each scheme, as done for steady-state 
flow, for example, by Leonard and Mokhtari (1990). A preliminary grid-refinement study 
of this kind is reported in Appendix 3, where it is shown that the uniformly third-order 
UTOPIA scheme is indeed significantly more efficient than the “quasi-third-order” scheme 
(which, in fact, is only second-order accurate). The stability region for the “quasi-third- 
order” method is the full square: (1 — |c Jt |)(l — |c y |) > 0. 

The advective-characteristic-tracking idea implied by (27) could, of course, be applied 
to uniformly higher-order interpolation methods. However, above third-order, algorithmic 
complexity increases rapidly, especially in terms of additional transverse and other mixed 
cross-difference terms. One strategy that appears to give quite good results without undue 
algorithmic complexity is to keep all third-order transverse and cross terms but expand the 
computational stencil for each advected face- value only in the normal direction. In other 
words, the one-dimensional contribution is expanded to higher order while keeping only 
third-order additional terms. For example, for the west face 

<t> w = <t> w ( higher -order ID) - °-l GRADTW 

- 5l[(l - c y )CURVTW + (1 - cJTWISTW] ( 59 ) 

This procedure has been explored by Leonard and Niknafs (1989), and appears to give good 
resolution and isotropy. Computational efficiency can be greatly enhanced by using an 
adaptive-stencil-expansion strategy; i.e., UTOPIA is used in smooth regions, whereas in 
regions containing large spatial gradients (or strong changes in gradient), the algorithm 
automatically extends the one-dimensional contribution to higher order on the basis of some 
non-smoothness monitor. This is highly cost-effective because the (more expensive) wider 
stencil is used only in isolated narrow (steep-gradient or high-curvature) regions involving 
a relatively small number of grid points. 

The main problem with multidimensional, collocated polynomial methods (above first 
order) is that they are not inherently positivity preserving. This is also true, of course, in 
the one-dimensional case. For a given continuous profile, higher-order methods are better 
in this respect, but for near-discontinuities, any polynomial scheme will generate spurious 
unphysical overshoots and undershoots, which may feed back and totally disrupt the solution. 
This phenomenon has been demonstrated for one-dimensional scalar advection (up to ninth- 
order) by Leonard (1991). In one dimension, simple flux-limiting strategies can be devised 
to guarantee monotonicity preservation. The two-dimensional extension of these ideas 
appears to be much more difficult. Some flux limiters (such as one-dimensional TVD 
conditions applied coordinate-wise) appear to be sufficient for maintaining positivity under 
some (rather restrictive) conditions, but tend to generate strong grid-dependent distortion. 
Some progress in devising truly multidimensional (and less shape-distorting) flux limiters 
has, however, been made. This is discussed in the next section. 
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MULTIDIMENSIONAL FLUX-LIMITERS 

Before tackling the multidimensional case, it is instructive to review the concepts of flux 
limiters for one-dimensional scalar advection. In this case, the model equation is 

it + u it = 0 (60) 

dt dx 

where u is a (positive) constant. An explicit single forward time-step, control-volume update 
formula for this equation can be written in compass-point notation as 

</>; = 4> P + c(4> w - <f> e ) ( 61 ) 

referring to west and east “effective” face-values (Leonard, 1991), where c = uAt/h is the 
Courant number. It is convenient to introduce a “generic” upwind-biased stencil. For 
example, near any given control-volume face (/), define upwind ( U ), downwind (D) and 
central (C) node values (on the basis of the sign of u) as shown in Figure 17(a). Now define 
a locally normalised variable 


4>(x) 


$(*) 4>u 

<$>d - <t>u 


(62) 


so that, in particular, <j> n v = 0 and <i> n D = 1 . Figure 17(b) shows the corresponding <f> n c node- 
value and the normalised face-value, 4> f . 

In one-dimensional flows, the main problem associated with maintaining monotonicity 
occurs when 4> p (= 4 > n c ) has a small (positive) value. If the interpolated 4> f value is too high, 
there will be too much flux advected out of the downwind face of the control-volume centred 
at C, and (= <^ +I ) may be drawn down below zero, thereby destroying monotonicity. 
A different interpolation could avoid this problem by predicting a smaller <j> f value. This 
is shown schematically in Figure 18, where a parabolic interpolation (a) for the face-value, 
4> f , is compared with an exponential interpolation (b), for the same 4> n c value. Over a time- 
step At, it is clear that the advective outflux predicted by case (a) would be larger than that 
of case (b). In order to derive the limiting conditions, assume first that 4> f lies between^" 
and 1 

4> n c < 4> f < 1 < 63 ) 

consistent with local interpolate monotonicity. Now consider (61) written in terms of 
normalised variables, and require the updated normalised value not to be negative 


K = 4> P + c(i w - 4> e ) > 0 (64) 

The situation is shown schematically in Figure 19. Referring to (63), written for the 
upstream (west) face, the worst-case (smallest) value for 4> w would be zero. From (64), this 
generates the condition 
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or, in general for any face, using the generic stencil of Figure 17, 

i f < — for 0 < i" < 1 ( 66 ) 

J c c 

One also needs to to consider suppression of overshoots when 4> p (= 4> c in Figure 18) 
is slightly less than 1 , for example. In this case, conditions may need to be imposed on the 
upstream-face influx to keep it small enough so that the updated 4>* p does not exceed 1 . But 
(63) is sufficient to avoid this because the largest ^-value is i p , and this is also the 
smallest ^ -value, so that 4 v can never overshoot its adjacent downstream value (at very 
worst, it remains constant over a particular time-step). Thus (63) and (66) taken together 
are sufficient to guarantee monotonicity -preservation (i.e., suppression of both overshoots 
and undershoots) in the case of one-dimensional advection at constant velocity; and this has 
been borne out by numerical experiments (Leonard, 1991). Figure 20 summarises these 
“universal” limiter constraints for the monotonic range, 0 < <j> n c < 1 . In the non- 
monotonic range, several strategies are possible, one of the simplest being to set 

4> f = 4> n c for 4> n c < 0 or 4> n c > 1 (67) 

this is also shown in Figure 20. The universal limiter is the basis for a number of higher- 
order, so-called “ULTIMATE”, one-dimensional advection schemes developed by Leonard 
(1991). The acronym stands for “universal limiter for transient interpolation modelling of 
the advective transport equations.” 

The flux-limiter constraints used by total-variation-diminishing (TVD) schemes are 
rather more restrictive than those shown in Figure 20. First, it should be noted that TVD 
schemes use second-order time averaging. Referring to Figure 21, this means that the face- 
value is assumed to vary linearly with time, so that 


4> f 




uAt 

~ 2 ~ 


= (1 -c)4> n f + c4> n c (68) 


recalling that 4>f (without a superscript) denotes the time-average over At , whereas 4>} is the 
interpolated face-value at time-level n. As explained by Leonard and Mokhtari (1990), the 
TVD flux-limiter constraints summarised, for example, by Sweby (1984) can be written in 
terms of <j> n f as 


if ^ 2 4> n c 

for 

0 

< 4> n c < 0.5 

(69) 

4> n f < 1 

for 

0.5 

< i n c < 1 

(70) 
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with 

0* > l n c for 0 < 4> n c < 1 ( 71 ) 

and 

tfy = 4> n c for j> n c < 0 or 0" > 1 ( 72 ) 

These are shown diagrammatically in Figure 22 (note that the vertical axis in this figure is 
4>", not 4> f ). In order to compare with Figure 20 (where the ordinate is 4> f , not 4y), it is 
necessary to rewrite the TVD constraints (69)-(72) in terms of 4> f , using (68). This gives 

4> f < (2 - c) 4> n c for 0 < 4> n c < 0.5 ( 73 ) 

and 

4> f < 1 + c(4> n c - 1) for 0.5 <j> n c < 1 ( 74 ) 

with the other constraints remaining the same. The corresponding diagram is shown in 
Figure 23, where the slopes of the upper constraints are given by 

5, = 2 - c (75) 

and 

S 2 = c (76) 

These should be compared with = 1/c and S 2 = 0 shown in Figure 20 for the (less- 
restrictive) universal limiter. The relationships between the upper-constraint slopes of 
ULTIMATE and TVD schemes are shown in Figure 24 as functions of Courant number, 
between c = 0 and 1 . 

In order to develop a two-dimensional positivity-preserving (PP) algorithm, one might 
try to use one-dimensional TVD or ULTIMATE fluxes in each coordinate direction sepa- 
rately (with c replaced by the local normal component Courant number at each face). For 
example, the TVD (CLAM) scheme of van Leer (1974) can be described in terms of <j> n f by 

4> n f = 2 4 > n c - (^ n c ) 2 for 0 < 4> n c < 1 ( 77 > 

with <j>" = <j> n c elsewhere. This is shown in Figure 25. Note how the slope of the parabolic 
curve is tangent to the TVD constraints as </>" -* 0 + and 1_. Also note that <£" = 0.75 when 
4>" c = 0.5, making this a second-order scheme in “smooth” regions (<j> n c -* 0.5). Figure 26 
shows the 2D results (for c x = c y = 0.25 and 9 = 45°) of using the ID CLAM scheme 
coordinate-wise. Although the scheme is evidently positivity preserving, there is con- 
siderable grid-orientation-dependent distortion, with cross-wind elongation and a tendency 
towards flattening along grid lines, producing a diamond-shaped distortion of the initially 
circular profile. The cross-wind elongation is primarily due to the lack of GRADT-terms. 
If these terms are included (before limiting), cross-wind elongation is reduced but the profile 
still tends to “square-up” along grid-lines. Also, as shown theoretically below, there are 
rather severe PP restrictions on the allowable component Courant-number ranges. 

Figure 27 shows results of using the one-dimensional third-order ULTIMATE 
QUICKEST scheme in each coordinate direction separately. In addition to cross-wind 
elongation (due primarily to the lack of cross-terms), it can be seen that this scheme is not 
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positivity preserving. There is no point in trying to improve isotropy because, as is shown 
theoretically below, using the one-dimensional universal limiter coordinate-wise cannot 
guarantee positivity. The following analysis is based on a suggestion by Bull (personal 
communication, 1990). Referring to Figure 28, consider a normalised monotonic “front” 
passing through a control volume so that post-frontal node-values are essentially zero, with 
pre-frontal values close to unity. At the time indicated, the positions of the leading-edge 
(LE) and trailing-edge (TE) of the front are such that the normalised node-value at the 
control-volume center is a small positive value, i p = l. This is the analogue of the one- 
dimensional case shown in Figure 18. If outfluxes significantly exceed influxes, the central 
node-value at the next time-level could be drawn down to negative values, thus destroying 
positivity. To avoid this, the following condition would be required for the updated 
variable, 4 > p , 

K = e + c x (&. w - i t ) + c y (4> s - 4> n ) > 0 ( 78 ) 

For the situation illustrated ( c x , c y > 0), the worst-case condition corresponds to <j> w and<£ f 
both being as small as possible (=0) and i e and <t> n both as large as possible. According 
to Figure 20 (with <j> n c = e), using the individual normal component Courant numbers, the 
universal-limiter constraints imply 


(large) = — (79) 

C x 

and 

(large) = 1 (80) 

c . 

Thus (78) becomes 

♦ c,(° - 1) * c,(. - 1) - -* < 8 » 

thereby violating the PP inequality. This means that the one-dimensional universal limiter 
used coordinate-wise cannot guarantee multidimensional positivity preservation. 

On the other hand, consider the TVD constraints of Figure 23. In this case (with 
4> w and 4> s still both equal to zero), 

^(large) = (2 - c x ) e ( 82 ) 

and 

4> n (large) = (2 - c y ) e ( 83 ) 

Hence, the updated central node-value now becomes 

fa - e + cJO - (2 - c x ) e] + c y [0 - (2 - c y ) e] 

= e [l -c x (2 -c x )-c y ( 2 — C y )] 


(84) 
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This expression is non-negative only if 

(c x ~ l) 2 + (c y - l) 2 > 1 for 0 < (c x , c y ) < 1 (85) 

This represents a rather severe restriction on the allowable ranges of the component Courant 
numbers, as seen by the circular arcs in Figure 29. In particular, for 9 = 45°, c x and c y 
must both be less than (1 - ^1/2) * 0.293 . Even more restrictive conditions are required 
in three dimensions (see Appendix 4). 


First Multidimensional Limiter 

In order to devise less restrictive positivity-preserving conditions, return attention again 
to Figure 28 and Equation (78). Assume that the slope 5! in Figure 20 is to be found in 
order to satisfy (78). Also assume that (63) still holds, so there can be no (normalised) 
overshoots (because, at least for a solenoidal velocity field, total influx cannot then exceed 
total outflux, for the same reasons as in the one-dimensional case) and, therefore, the only 
condition that needs to be considered is the avoidance of undershoots. This means that, 
instead of (79) and (80), the worst-case outflow face- values are written in terms of the (as 
yet unknown) slopes, S? and S y , as 

^(large) = S? e ( 86 ) 

and 

<£ n (large) = S? e ( 87 ) 

where (for the moment) it is assumed that Sf and S? may be different but related by 
symmetry according to 

s; * s,(c,, «,) < 88 > 

and 

S? = Sfc y , c x ) ( 89 ) 

Substitution into (78) gives 

K = 6 + c x (0 - S? e) * c y (0 - S? e) > 0 ( 90 ) 

resulting in the condition 

1 - c Si - c Si y > 0 ( 91 ) 

** / 

The most general solution of this inequality would require 


and 


Si < 


f( c *> c y ) 


c x f(c x , c y ) + c y f(c, c x ) 


= S l( C *> C y) 


(92) 


Si < Sfc , c x ) 


(93) 


where / is a (not necessarily symmetrical) function of c x and c y , to be chosen. Note that if 
c y = 0 and/(0, c x ) is finite, then 
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S,(c,, 0) = 1 (94) 

c, 

consistent with the one-dimensional formula. Also, at 0 = 45°, 

S,(c, c) = -L (95) 

2c 

This means, in particular, that at c x = c y = 0.5, 

S,(0.5, 0.5) = 1 (96) 

regardless of the form of/. If/ is in fact symmetrical, then it cancels in (92) leaving 

S,(c„ C y ) = S? = S, * - 1 - - (97) 

k,l * k,l 

where absolute values have been included to generalise to cases involving negative velocity 
components. This is the least restrictive (i.e., largest) symmetric form satisfying (91). Less 
restrictive non-symmetrical forms might be found by specifying the functional form of /. 
However, this strategy will not be pursued here. Of course, S, cannot be less than one, as 
this would violate the left-hand inequality of (63). Figure 30 shows contour lines of S { 
within the allowable region (5! > 1). Results of using this two-dimensional limiter in 
conjunction with the QUICKEST scheme in each coordinate direction are shown in Figure 
31 for 0 = 45°. The simulation is positivity preserving, but again there is considerable 
cross-wind elongation due to lack of cross-terms, together with steepening in the flow 
direction and concomitant flattening of the peak. 

One could, of course, restore the cross-terms missing from the QUICKEST scheme and 
apply the limiter, (97), to the fully-third order UTOPIA scheme. Results are shown in 
Figure 32 for c x = c y = 0.25. The diamond-shaped distortion is quite disappointing; as 
advection proceeds further, the profile tends to square-up along grid lines. Apparently, the 
limiter represented by (97) is overly restrictive. For example, (97) shows that the slope of 
the limiter becomes 1 along the boundary lines 

\c,\ ♦ |c,| = 1 (98) 

seen in Figure 30. This means that, for Courant numbers in ranges near these boundaries, 
the simulation is approaching the highly artificially diffusive and anisotropic first-order 
scheme given by (19), with results similar to those shown in Figure 4. 
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Second Multidimensional Limiter 

Some progress toward better isotropy can be made by incorporating the transverse- 
gradient terms after the limiting step. Assume that is a higher-order face-value not 
including GRADT terms. For example, the west-face UTOPIA 4> f would be given by (58), 
omitting the GRADTW term; i.e., 

o c z (1 - c]) 

<l> w = <t> UN - — GRADNW - — CURVNW 

2 6 

- ^[(1 - c y ) CURVTW + (1 - c x ) TWISTW] 

Now assume that this temporary face value, written in normalised form as 
constrained, near 4> n c -* 0 t , by 

4> n c < 4>° f < 0 REF1 = 5 , 4> n c 

where remains to be found in order to satisfy PP conditions. Working in terms of 
w/mormalised variables, the right-hand reference value can be written 

4>REF1 = <*>v + ( 4> c - K) < 101 ) 

using the generic-stencil notation of Figure 17. 

Referring to Figure 11, the update algorithm in w/mormalised variables can be written 
(for constant velocity) as 


(99) 
0°, is 

(100) 


4>p = <t> P + c x (<t> w - <f > e ) + c (4> s - <t > n ) 


( 102 ) 


Now write the actual face-values as the algebraic sum of the temporary values and the 
appropriate transverse-gradient terms. For c x and c y both positive, these take the form 


K = 4>l ~ -f (K - 


(103) 


' *s) 


(104) 


*S - - y (<t>s - *sw) 


(105) 


<t>n ~ 4>n ~ y (4> P ~ 4> w ) 


Substitution into (102) gives 


<t>p ~ <t> P + c x c (4> P - <t> w - (t> s + <t> sw ) + c x (<t>° w - 


<t>°<) * c y (<t>° s 


* 2 ) 


(106) 


(107) 
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Consider first, non-negative node- values increasing in the flow direction, with 4> p quite 
small. The aim is to avoid the updated <t>* P becoming negative. Once again consider worst- 
case small influxes and large outfluxes. If <fr w and <j> s are small (but non-negative), then the 
left-hand inequality of (100) gives the worst-case inflow values as 

<£°(small) = <f> w (108) 

and 

$°(small) = <£ 5 (109) 

The largest outflow values need to be constrained by the slope S x . In terms of unnormalised 
variables, referring to (101), 


<*>°(large) 

= <t>w + (<t>p ~ ^w) 

(110) 

4>°(large) 

= 4>s + s l ( <t> P - <t>s ) 

(HI) 


where again (for the moment) S l and Si are assumed to be possibly different but symmet- 
rically related. Substitution of (108)-(1 1 1) into (107) gives 

4>p = <t> P [l + - c x Si - c y 5/] + 4> sw 


+ c x (Si - c y ) 4> w + c y (Si - c x ) <t> s 


( 112 ) 


and this is required to be non-negative for arbitrary, non-negative 4> P , <f> sw , (f> w , and <t> s . 
Thus the coefficients of each of these node-values must all be non-negative. An upper bound 
on Si and Si is obtained by setting the coefficient of <£ P to zero, the corresponding 
conditions for <f> w and <t> s providing lower bounds. The upper bound is chosen, as this yields 
the least restrictive limiter. This implies 


S 


x 

1 


< 


and, by symmetry, 


(1 + c x c y )f(c x , c y ) 
c ,/( c *> c y ) * c y f(c y , c x ) 

Si < S x (c y , c x ) 


Si(c x , c y ) 


(113) 

(114) 


where, again, / is some (in general not necessarily symmetrical) function to be specified. 
As before, this 5, is consistent with the one-dimensional case; but at 45° 


5,(c, c ) 


1 + c 2 
2c 


(115) 


which is now always greater than (or equal to) 1 . 

Once again, if/ is a symmetric function of c x and c y , the (only) symmetrical form of 
Si becomes 
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Si = S? < 5, = 

From (112) this gives 

4>; > c,c y <t> sw 


(! + C x C y) 


( C x + o 


^(1 - O 


<t>w 


cM ~ c l x ) 


4 ><. 


(116) 


(117) 


( C , + C y) (S + O 

And since, for stability, \c x \ and |Cy| are both less than or equal to 1, the updated node- 
value is indeed non-negative. Unsymmetrical forms could perhaps lead to somewhat less 
restrictive conditions, but these are not pursued here. 

However, it turns out that (116) is not sufficient for giving PP results in general. One 
must also consider (positive) node-values decreasing in the flow direction, as well. This will 
require restrictions on 4> f near 4> n c -* 1.. Specifically, assume that, in normalised-variable 
notation, 

(118) 


70 


near <t > " 


<t>c * <*>; S ^REP 2 - 1 + S 2 (K ~ 1) 

1_. In terms of «/mormalised generic variables, (118) can be rewritten as 


<t>c > 4>r 


4> 


REF2 


= + S 2 ( 4> c - <i> D ) 


(119) 


using the notation of Figure 17. The inequality reversal is due to the negative normalising 
denominator in this case. In a decreasing region, worst-case conditions (tending to make the 
control-volume central node-value negative) imply 


4>e (large) = <t> F 


and 


</>„( large) = 4> p 

as well. The smallest allowable inflow face-values, according to (119), are 

0® (small) = 4> P + S 2 0 <t> w - 4 > P ) 


and 


where 


and, by symmetry, 


<t>°s( small) = <t> p + S 2 ( 4> s - <t> P ) 


S 2 = S l( c x > c y) 


S 2 y = S 2 (c , cj 


( 120 ) 

( 121 ) 

( 122 ) 

(123) 

(124) 

(125) 


Substitution of (120)-(123) into (107) gives 
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<t>; = Ml + C x C y - C x S 2 ~ C y S 2) + C x C y M 

+ C x 0^2 - C y) 4> w + Cy ($2 - c x )4>s 

Choosing the smallest (and therefore least restrictive) S 2 ’s as 


and 


S, = c 


Si = c x 


eliminates the <j> w and 4> s terms and results in 


(126) 

(127) 

(128) 


<t>; = Ml - C x c y) + C x C y <i>SW 


(129) 


Since this represents a weighted linear interpolation between the two (nonnegative) node- 
values, <l > sw and <f) P , the result (for \c x c y \ < 1) is also nonnegative, approaching <^> F for small 
\c x c y \ and equalling <f> sw for \c x c y \ = 1 (consistent with exact point-to-point transfer across 
the grid diagonal). 

The two-dimensional PP limiter restrictions on are summarised in Figure 33, where, 
for a constant velocity and corresponding component Courant numbers, c x and c y , allowing 
for positive or negative velocities, 


and 


S t 


1 + \c x c y \ 

kj + \Cy\ 


So = 




(130) 

(131) 


Contour plots of (130) are shown in Figure 34. As seen from the contours in Figure 34, 
there is a rather large region where 5, is only slightly greater than 1 (e.g., see the large 
region between the contours 5, = 1.25 and 5, = 1); however, this is still less restrictive than 
the first two-dimensional limiter, given by (97). 

Figures 35 through 39 show results of using the limiter given by (130) and (131) in 
combination with the UTOPIA scheme (i.e., with transverse-gradient terms restored after 
limiting but all other terms included beforehand). In Figure 35, c x = c y = 0.25 (0 = 45°); 
in Figure 36, c x = 0.25 and c y = c x /2 (0 * 26.565°); in Figure 37, c x = 0.25 and 
c y = c x /3 (0 ® 18.435°); and in Figure 38, c x = 0.5, c y = 0 (0 = 0°). These results are 
seen to be generally much better than those using the first two-dimensional limiter given by 
(97). However, at c x = c y = 0.5 (0 = 45°) shown in Figure 39, there appears to be some- 
what more distortion. In this case, the limiter slopes are S j = 1.25 and S 2 = 0.5, repre- 
senting rather restrictive constraints. At larger Courant number values the simulation 
becomes noticeably more diffusive. 
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Extension to Variable Velocity 

The two-dimensional flux-limiters developed so far have been based on a uniform 
advecting velocity field at an angle to the grid. Clearly, this is too restrictive to be of use 
in practical calculations involving variable velocity fields. However, the same principles can 
be generalised to variable solenoidal velocity fields satisfying the constant-density continuity 
equation 

V-V =0 ( 132 ) 


For simplicity, only the first two-dimensional limiter will be extended here. The extension 
of the second limiter is more complicated and will be taken up in a future report. 

Consider the update equation for variable velocity using conventional compass-point 
notation at each face 


<b* P - <f>„ + c <t> -c6+C(t>-c<t> 

^ r 'r p '■'xw xe ^ e ys ^ s yn 


(133) 


where, by continuity 




c + 

xe 


c - c 

ys yn 


= 0 


(134) 


Four separate cases need to be considered: 

(i) Two inflows on adjacent faces together with outflows on the other faces. This 
is represented by all c’s being positive, for example. 

(ii) One inflow and three outflows. 

(iii) Three inflows and one outflow. 

(iv) Inflows on opposite faces and outflows on the remaining faces. 

These are sketched in Figure 40. In Case (i), Figure 40(a), the worst-case update becomes, 
in terms of unnormalised variables, 


<f>* P = <f> p + c xw 4> w - c xe [<t> w + S, (4> p - <f> w )] + c ys 4> s - c yn [<i> s + S { (<f> p - 4> s )] ( 135 ) 

or 

= *,[i - s,(c„ * c„)] - kK + - o] ♦ * c >- W - *>] < 136 > 

Choosing 

S . = (137) 

+ Cy. 

gives positivity-preserving results. This is the appropriate limiter slope for use on both the 
east and north faces— i.e. , on the outflow faces for the CV cell under consideration. For the 
conditions of Case (i), but for other combinations of signs of the Courant number 
components, (137) can be immediately generalised, for two outflowing faces, to 




1 



(138) 


using an obvious notation. 
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For Case (ii) shown in Figure 40(b), the worst-case update becomes 
<t>p = <t>p + c m 4> w - C xe [<t> w + 5, ( <t> P - <t> w )] 

~ \ C ys\ [^N + S l (4>P ~ ^v)] ~ C >»[^5 + S l (<t>P ~ ^s)] (139) 

K = <*> p [i - S x {c xe + + \c ys \)] * kK + C ~( S 1 - D] 

+ - 1)] + <t> N [\c ys \ (5, - 1)] (140) 


In this case, noting that c xe and c yH are both positive, 




1 


+ c. 


yn 




(141) 


is the appropriate limiter slope for each of the three outflow face-values. Again, this can be 
generalised for three outflowing faces as 


Sx 



(142) 


It is not difficult to see that similar formulae would apply in Cases (iii) and (iv). 

In summary, the limiter slope for any face is computed as follows: 

(i) Identify the control-volume cell for which the face in question is an outflow 
face. 

(ii) Take the S r slope for that face as 




£ k 


out I 
n 


(143) 


where the sum is over all outflow faces for the identified CV cell (including the 
face in question). 

A similar strategy can easily be devised for three-dimensional flow in an obvious way. 
This type of flux-limiter has been tested for scalar advection in 2D and 3D solenoidal 
velocity fields. As expected, results are positivity-preserving to machine accuracy, even 
when discontinuities are present in the initial conditions. Details will be given in a 
subsequent report. 
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CONCLUSION 

Advection is a key dynamical process in many areas of fluid dynamics. As computer power 
has increased and the spatial resolution of CFD codes has improved, so the need for more 
accurate numerical advection schemes has increased in significance. This is because the 
higher spatial resolution allows (in principle) the representation of sharper gradients and 
changes in gradient. More detailed resolution of advective transport means that less 
empiricism has to be used in modelling important fluid-dynamic processes. Unfortunately, 
as is well known, traditional advection methods, such as the space-time-centred leap-frog 
scheme, do not perform well under highly advective conditions, often exhibiting large 
regions of spurious numerical oscillations. Such behaviour may be manifested in CFD 
simulations as, for example, negative values of density concentrations, or other unphysical 
phenomena. Numerically dispersive advection schemes of this type may even generate a 
kind of “numerical turbulence” that, in turn, can lead to the erroneous prediction of 
anomalous mixing processes. These features may interact with the physical parameter- 
ization within the models and lead to significant degradation of numerical simulations and 
forecasts. 

Two aspects of accurate advective modelling have been addressed in the present report. 
The first concerns the problem of anisotropic distortion— i.e., the dependence of the solution 
on grid orientation. Methods based on using one-dimensional formulae in each coordinate 
direction separately show strong anisotropic distortion at non-zero flow-to-grid angles. The 
introduction of transverse-gradient cross-terms into the advective fluxes greatly alleviates this 
problem in second-order methods. However, the multidimensional, uniformly second-order 
method (the generalised conservative Lax-Wendroff scheme) is still too numerically 
dispersive for use as the basis for a good CFD advection code. Fortunately, the techniques 
used to derive this method generalise to higher order. As in the case of one-dimensional 
advection, the multidimensional, uniformly third-order, polynomial interpolation algorithm 
(UTOPIA) possesses several desirable properties. The conservative, control-volume 
formulation is relatively straight-forward algorithmically, isotropy appears to be excellent, 
and phase error is low. Because of the flux-based, conservative formulation, “mass” is 
conserved (to machine accuracy). The only drawback (as with all collocated polynomial 
schemes above first order) is that it is not inherently positivity-preserving. 

In order to produce a multidimensional, positivity-preserving advection scheme, one 
might consider using well-known one-dimensional TVD flux-limiters in each coordinate 
direction. However, this approach has been shown here to have two shortcomings: the 
allowable Courant-number range is very limited, and results are again highly grid-dependent, 
showing a strong “squaring-off” distortion along grid-lines. The anisotropy appears to be 
due partly to the use of one-dimensional limiter formulae (which are inherently grid- 
orientation dependent) and partly to the fact that TVD limiters, in general, are overly 
restrictive, tending to add (coordinate-wise) artificial diffusion locally, in order to suppress 
potential oscillations. 

One of the main focuses of the current work has been to generalise the less-restrictive 
universal limiter concepts originally developed for one-dimensional flows by Leonard (1991). 
The one-dimensional universal limiter applied coordinate-wise in two dimensions has been 
shown theoretically not to preserve positivity. However, the analysis used in deriving this 
result also suggests two different multidimensional generalisations that are symmetric with 
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respect to x- and y-component Courant numbers (and are consistent with the one-dimensional 
case). Unfortunately, the first multidimensional limiter, (97), reintroduces considerable 
anisotropic distortion. The second multidimensional limiter, given by (130) and (131), is 
based on adding the transverse-gradient terms after the limiting process. In this way, 
reasonably good isotropy is maintained when applied to an inherently isotropic advection 
scheme such as UTOPIA. Fortunately, both the advection scheme itself and the limiter con- 
straints can be generalised in a straight-forward manner to three dimensions (see Appendix 
1 and Appendix 5, respectively). Appendix 2 also outlines the inclusion of physical diffusion 
terms to a consistent order, and a simple upwinding strategy designed to avoid IF statements 
is sketched in Appendix 6. 

The multidimensional third-order advection scheme and the generalised limiter 
constraints were initially derived on the basis of uni-directional constant-velocity linear 
scalar advection. Similar principles have been applied in the case of spatially varying 
solenoidal velocity fields using the first multidimensional limiter. Generalisation of the 
second multidimensional limiter to variable velocity fields is more complicated, and will be 
discussed elsewhere. The main aim of this research project has been to advance the state 
of the art of advection modelling beyond standard techniques (generally involving gross 
anisotropic distortion and potentially disastrous numerical dispersion). And, to the extent 
that this has been achieved, highly advective modelling should be enhanced. This will 
become apparent when the current strategy has been implemented in production codes and 
applied to real fluid-dynamic processes. 
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APPENDIX 1 

THREE-DIMENSIONAL FORMULAE 

In three dimensions, the computational molecule for third-order upwinding consists of the 
twenty nodes shown in Figure A. 1.1, where the notation refers to the conventional compass- 
points ( N , S, E, W) together with top ( T ) and bottom ( B ) levels. Note the upwind bias, 
assuming that the Courant number components, c x , c y , and c z are all positive. Taking local 
coordinates centered at P, the </>-value at time-level n can be written 

v(z,v,n = c, + c 2 £ + c 3 e + c,e 

+ C 5 1) + C 6 V 2 + C 7 T) 3 

+ c 8 r + c 9 r 2 + c 10 r 3 

+ Cntv + C l2 (v 2 * C 13 Z 2 V 

+ J7 f + c 15 r) i 2 + c 1# *? 2 f 

+ c 17 r? + c 18 ri 2 + c 19 r 2 ^ + c 20 ^r (a.i.d 

This should be compared with (44); note the obvious similarity in structure, but also note 
the appearance of the triple-product C 20 -term. Because of the similarity with the two- 
dimensional case, evaluation of the coefficients by collocation follows a predictable pattern: 

C, = 4> p (A. 1.2) 


~ K _ <ft £ ~ 3 <t>p + 3 <f> w - 
2 Ax 6 Ax 


(A. 1.3) 


C 3 


<t>E ~ 2 4> P + K 

2 Ax 2 


(A. 1.4) 




- 3<fr p + 3<fr w - 
6Ax 3 


(A. 1.5) 


~ 4>S _ ~ 3< ^P + 3 4>S ~ ^55 

2 Ay 6 Ay 


(A. 1.6) 


C 6 


<t> N - 2 4> p + <t> s 
2 Ay 2 


(A. 1.7) 


30 


C 7 = 


<t> N - 3<t> P + 3 <t> s - <t> ss 


6Ay 2 


C. = 


&T ~ &B _ &T ~ 3<t> p + 3 <t> B ~ <f> 


BB 


2Az 


6Az 


C 9 = 


<t> T ~ 2<t> p + <t> B 

2 Az 2 


'u 


'14 


'17 



C in 

4> t 

- 3 <t> p + 3 <t> B - 

^BB 




10 


6 Az 3 




<** 

-K) 

~ Wse 

- 0 S ) + (<£„ - 

^Nw) 

- (4>p - 

K) 




2 Ax Ay 




' = 

(4>„ - 

■ 4>nw) 

- 2(<t> P - K) 

+ (*5 

- 


12 



7 Ax Ay 2 




’ = 

(4> e - 

■4>se) 

- 2(*, - *,) 

+ (<*V 



13 



2 Ax 2 Ay 




(<t>E 

- 4> P ) 

- (<t>BE 

_ <£/») + (</>7- " 

^7w) 

- (<£/> _ 

<t>w) 




2 Ay Az 




\ c = 

(<t>T ~ 

4* TW ) 

- 2(*, - *») 

+ (4> S 



15 



2 Ay Az 2 




= 


‘ 4> BN ) 

~ 2 (4> p - (f > B ) 

+ (4>S 

- *«) 


'16 



2 Ay 2 Az 




(4>£ 

-K) 

- (<*«S 

~ 4> b ) + (<£ r " 

^7w) 

~ (<£/> _ 

K) 


2AzAx 


r _(</>£- <I>BE ) " 2 ^P - 4> B ) + “ ^Bw) 

'-'I 


"18 


(A. 1.8) 
(A. 1.9) 
(A. 1.10) 
(A. 1.11) 
(A. 1.12) 
(A.1.13) 
(A. 1.14) 
(A. 1.15) 
(A. 1.16) 
(A. 1.17) 
(A. 1.18) 


2 Az Ax 2 


(A. 1.19) 
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(<t>T ~ ^7Tv) ~ 2(</>p ~ <ft w ) + (<fr B ~ <ft flw ) 

2 A z 2 Ax 


(A. 1.20) 


(4>f» - 4>,y - <ft 5 •*• ^5tv) ~ ~ Kv ~ + ^«sw) (A. 1.21) 

Ax Ay Az 


In terms of a three-dimensional velocity, V = ( u , v, w), the advective back-tracking 
formula can be written 


When this and the above C-values are substituted into (A. 1. 1), the resulting formula can be 
rearranged into a flux-based, conservative, control-volume, single time-step, explicit update 
algorithm of the form 

4>* P = + FLUXW(i, j, k) + FLUXS (/, j, k) + FLUXB(/, j, k ) 

- FLUXW(/'+l, j, k) - FLUXS (/, ;+ 1, k) - FLUXB(/, j, k+1) (A. 1.23) 
where, for example, 


and so on. By analogy with the two-dimensional case it is not difficult to show that the 
advected west face-value, for example, is given by 


<f>p = (f> n (-uAt, -vA t, -wAt) 


(A. 1.22) 


FLUXW(t, j, k ) = c x <t> w 


(A. 1.24) 


c x (i - c;) 

<£ w = <i> UN - GRADNW - — CURVNW 

2 6 


GRDTYW - — GRDTZW 


2 


2 


C -1 [(1 - c y ) CRVTYW + (1 - c x ) TWSTYW] 
C -± [(1 - c z )C RVTZW + (1 - c x ) TWSTZW] 



(A. 1.25) 


Where <1> UN is the linear interpolation across the west face 



(A. 1.26) 
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the normal gradient is given by 


GRADNW = 4> P ~ 4> w 

(A. 1.27) 

and the upwind-biased normal curvature is 


CURVNW = 4> P - 2 4> w + 4>w 

(A. 1.28) 

just as in the one-dimensional formula. The two (two-dimensional) upwind-biased trans- 
verse-gradient terms are given by 

GRDTYW = 4> w ~ <t> sw 

(A. 1.29) 

and 

GRDTZW = <*> w - cf> BW 

(A. 1.30) 

The transverse-curvature and (x,y)- and (x,z)-twist terms also have the same (upwind-biased) 
form as their two-dimensional counterparts 

CRVTYW = </> w - 24> w * <t> sw 

(A. 1.31) 

CRVTZW = - 24> w + <f> BW 

(A. 1.32) 

TWSTYW = 4> p - 4> s - <t> w + 4> sw 

(A. 1.33) 

and 

TWSTZW = 4> p - <f> B - 4> w + <{> BW 

(A. 1.34) 


Finally, the last term in (A. 1.25) is a peculiarly three-dimensional term representing an 
upwind-biased (y,z)-twist contribution 

TWSYZW = -*„-*„* * m (A. 1.35) 

Figure A. 1.2 shows the ten node-values involved in computing 4> w (for c x , c y , c z > 0). The 
other face-values, 4> s and </> fc , can be obtained by using a straight-forward permutation 
process. However, this is not necessary, because a “generic” stencil can be used for any 
control-volume face, determined by the signs of the individual local Courant number 
components. This is outlined for the two-dimensional case in Appendix 6. The three- 
dimensional extension is straight-forward. 
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APPENDIX 2 

INCLUSION OF DIFFUSION TERMS 


Consider the ad vection -diffusion equation written in conservation form as 

0* + V-(Vtf>) = V • (DV<t>) 


(A.2.1) 


A numerical update formula can be written in the form of (A. 1.23) provided each flux con- 
tains both advective and diffusive contributions. For example, the west-face flux is given 
by 


FLUXW <«,/,*) - c„*. (A.2.2) 


introducing the diffusion parameter, a xw , defined by 


or 


Ax 2 


(A.2.3) 


Similarly for the other faces. The task now is to estimate the face-value and face-gradient 
appearing in (A.2.2). At third order (and above) one needs to consider the effect of diffu- 
sion in estimating <£ w and, conversely, the effect of advection in estimating (84>/dx) w . 

In order to assess these effects, consider the constant-coefficient one-dimensional 
advection-diffusion model equation 

= - u ?± + (A. 2. 4) 

dt dx dx 2 


Now formally integrate this equation in time, from 0 to r. At the west face, for example, 
this gives 





(A.2.5) 


The first integral on the right represents the variation of 4> w with time due to purely advective 
effects (obtained by setting D = 0); the second integral represents the diffusive contribution. 
The latter term is estimated simply as 


(>£*) dr = 

flD^fdr . 

Id r 

(A. 2. 6) 

' dx 2 ' w 

■'o ' dx 2 '* 

' dx 2 !» 



Averaging <f> w (r ) over r = 0 -*■ At, then leads to 


<*> 


w 



w 



<M adv) 


w 


At 

T 


(A.2.7) 
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where the first term is the purely advective contribution to the time-average and the second 
term an estimate of the diffusive contribution. Equation (A. 2. 7) can be rewritten as 


K = *„( adv) + ^ CURV B W 


(A.2.8) 


where CURV" is an appropriate second-difference. By equating CURV" with the upwind- 
biased (normal) curvature term, the one-dimensional third-order time-averaged face-value 
can be written 


= <*> 


UN 


- — GRADL - 
2 


(1 - cl) 


a 


X w 

T. 


CURVL 


(A.2.9) 


Clearly, (A. 2. 7) can be generalised to two and three dimensions by writing for the west 
face, say, 


<f> w = 0 w (adv) + (DV 2 <t>) 


n At 7 


(A.2.10) 


This generalises (A. 1.25) to the following form 

2 


<t>W = 


UN 


— GRADNW - 
2 


(1 - O _ 

6 2 


CURVNW 


fn GRDTYW - GRDTZW 
2 2 


c a 

jun - c ) - — 
4 yw 2 J 


CRVTYW - 1 - c x J TWSTYW 


c a 

— (1 - c ) - — — 
4 iw 2 


CRVTZW - _|T (1 - C„) TWSTZW 


CywCzw TWSYZW 


(A. 2. 11) 


For estimating the time-averaged normal gradient across a given face, refer first to the 
one-dimensional case sketched in Figure A.2.1. In time Ar, the curved profile shown is 
advected to the right; thus, the face-gradient is continuously changing with time according 
to the linear relationship 



(A.2.12) 


The time-averaged gradient over At is therefore 

(W = / a^\° _ /avr + 

\dxl \dxl \dx 2 / 


w 


2 


0(At 2 ) 


(A.2.13) 
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This can be rewritten as 


(i^l Ax * GRADW - fncURV^ (A. 2. 14) 

' dx I w 2 

where CURV^, is an appropriate second-difference in the vicinity of the west face. Once 
again it is convenient to use the upwind-biased curvature term, CURVW. In the one- 
dimensional case the approximate gradient is therefore given by 

(^) Ax = GRADW - — CURVW (A. 2. 15) 

\ dx lw 2 

In two and three dimensions, the normal gradient will vary with time due to twist terms, as 
well. It is not hard to generalise (A. 2. 15) for the west face, for example, to give 

l^] Ax = GRADNW - 5^ CURVNW - Sn TWSTYW - ^TWSTZW (A. 2. 16) 
\dxL 2 2 2 

The total west-face flux is given by substituting (A. 2. 11) and (A. 2. 16) into (A. 2. 2). 

For a model constant-coefficient form of the governing advection-diffusion equation 

0$. + \>V<t> = DV 2 <£ (A. 2. 17) 


it is possible to derive an exact three-dimensional complex amplitude ratio, G exict . Use of 
(A. 2. 11) and (A. 2. 16) and corresponding formulas for the other faces then produces an 
algorithm whose G- expansion agrees with that of G ex , ct through all third-order terms in the 
wave-number components. Specifically, consider a three-dimensional plane wave of the 
form 

4> = A(0 exp(i k*x) (A. 2. 18) 

where the vector wave-number is given by 

k = kj i + k y j + k z k (A. 2. 19) 

Substitution into (A. 2. 17) then gives 

<t> = /1(0) exp[-k*kDr] exp[t k*(x-Vr)] (A. 2.20) 


Assume constant grid spacing, Ax = Ay = Az = h, and define 

D At 


(A. 2.21) 


0 


kh 


(A. 2. 22) 
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and 


c = 


\At 


(A.2.23) 


Then the exact three-dimensional complex amplitude ratio is given by 


cxict 


= exp(-a 0 2 ) exp[-t(c*0)] 


where 0 2 = 0-0. Expanding real and imaginary parts gives 


G c«ct - 1 " 


a 0 2 + - (c • 0) 2 
2 


(A. 2. 24) 


+ %0 2 (cO) 2 + ±(c-0) 4 
2 2 24 


- i 


c-9 - 


ad 2 (c-fi) + I (c • 0) 3 
6 


— 0*(c-0) + “0 2 (c*0) 3 + _L(c-0) 5 
2 ' 6 v 120 


(A. 2. 25) 


A von Neumann analysis of the generalised UTOPIA scheme, incorporating the effects of 
diffusion on advection and vice versa, as described above, leads to 


G UTOPIA 1 


a 0 2 + - (c • 0) 2 
2 


- i 


c-0 - 


a0 2 (c-0) + i(c-0) 3 

6 


(A.2.26) 


thus matching G ex>ct through all third-order terms in the ^-components. 
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APPENDIX 3 

GRID-REFINEMENT STUDIES 

Some preliminary grid-refinement studies have been made for some of the schemes discussed 
in this report. All of the tests are for an angle of 0 = tan' ] (l/2) * 26.525°, using Courant 
numbers c x = 0.5 and c y = 0.25. The number of time-steps, N, is chosen so that the exact 
solution would move 2 units in the x-direction and 1 unit in they-direction. This means that 

N = — (A.3.1) 

c * h 

In addition to the base case of h~ l = 31, the grid is refined to h~ l = 61, 121, and 241. 

At each grid point, the final-time error is computed by subtracting the theoretically 
exact value from the computed value 

= hi ~ < M*,< yj) 

An Lj global error diagnostic is also computed from 

1/A 1/A 

i, - a'EEkJ 

.=.1 j-i 

The following tables show the total number of time-steps, N, the grid-refinement, h'\ 
the minimum nodal 0-value, the maximum nodal 0-value, the maximum (signed) error, the 
L x error, the CPU time (on a CRAY-YMP, to the nearest second), and the computation time 
relative to the 2D Lax-Wendroff scheme. The convergence rates for the maximum error and 
the L x error are taken to the nearest integer. 


(A.3.2) 
(A. 3. 3) 


Table A.3.1 2D Lax-Wendroff 


N 

} i' 1 

^ min 

4 > max 

^max 

u 

CPU time 

Rel. time 

124 

31 

-1.13 x KT 1 

0.822 

3.75 x 10' 1 

3.09 x 10' 2 

3 

1 

244 

61 

-1.43 x 10- 2 

0.965 

1.46 X 10' 1 

8.34 x lO' 3 

24 

1 

484 

121 

-3.16 x KT 6 

0.995 

3.85 x 10' 2 

2.13 x 10 43 

188 

1 

964 

241 

+2.16 x KT 11 

0.9994 

9.53 x 10' 3 

5.35 x lO* 4 

1485 

1 

CONVERGENCE RATE 

0(h 2 ) 

0(h 2 ) 
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Table A.3.2 ID QUICKEST with GRADT terms 


N 

h' 1 

$min 

^max 

^max 

u 

CPU time 

Rel. time 

124 

31 

-3.24 x 10- 2 

0.844 

-1.66 x lO' 1 

1.26 x lO' 2 

4 

1.208 

244 

61 

-2.19 X 10- 4 

0.973 

-3.91 x 10' 2 

3.20 x 10' 3 

29 

1.208 

484 

121 

-5.17 X 10‘ 8 

0.997 

-8.50 x 10' 3 

7.96 x 10' 4 

227 

1.208 

964 

241 

-1.43 x IQ' 11 

0.9996 

-2.02 x 10' 3 

1.99 x lO' 4 

1794 

1.208 

CONVERGENCE RATE 

0(h 2 ) 

0(h 2 ) 



Table A.3.3 UTOPIA 


N 

h' 1 

^min 

^max 

^max 

U 

CPU time 

Rel. time 

124 

31 

-5.93 x 10' 3 

0.872 

-1.28 x 10' 1 

6.47 x 10' 3 

4 

1.333 

244 

61 

-1.02 x 10' 5 

0.977 

-2.29 x 10' 2 

9.95 x 10‘ 4 

31 

1.333 

484 

121 

+6.01 x 10' 11 

0.997 

-3.07 x 10' 3 

1.30 x 10' 4 

249 

1.333 

964 

241 

+ 1.67 x 10' 11 

0.9996 

-3.86 x lO' 4 

1.60 x 10' 5 

1980 

1.333 

CONVERGENCE RATE 

0(h 3 ) 

0(h 3 ) 



Table A.3.4 UTOPIA with first 2D limiter 


N 

r 1 

& min 

<£m.x 

^max 

u 

CPU time 

Rel. time 

124 

31 

7.80 X 10' 8 

0.716 

3.95 x 10' 1 

3.15 x 10' 2 

5 

1.505 

244 

61 

2.48 x 10‘ 7 

0.888 

4.23 x 10' 1 

2.14 x lO' 2 

36 

1.505 

484 

121 

3.76 x 10' 8 

0.943 

2.95 x 10' 1 

1.31 x 10' 2 

283 

1.505 

964 

241 

1.75 X 10' 9 

0.956 

1.65 x 10’ 1 

7.40 x 10' 3 

2235 

1.505 

CONVERGENCE RATE 

0(h) 

0(h) 



Table A. 3.5 UTOPIA with second 2D limiter 


N 

K x 

& min 

4>m«x 

^max 

L\ 

CPU time 

Rel. time 

124 

31 

3.00 x 10' 9 

0.711 

-2.89 x 10' 1 

1.05 x 10* 2 

5 

1.623 

244 

61 

1.44 x lO' 10 

0.886 

-1.15 x 10' 1 

3.27 x 10' 3 

39 

1.623 

484 

121 

5.42 X lO' 11 

0.957 

+5.72 x 10' 2 

8.26 x 10~ 4 

305 

1.623 

964 

241 

3.35 x lO' 11 

0.984 

+2.59 x 10' 2 

2.06 x 10' 4 

2410 

1.623 

CONVERGENCE RATE 

0(h) 

0(h 2 ) 
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APPENDIX 4 

THREE DIMENSIONAL TVD CONSTRAINTS 

In three dimensions, TVD flux-limiters can be used component-wise but, as in the two- 
dimensional case, there will be restrictions on the component Courant numbers mandated by 
positivity preservation. The three-dimensional update equation corresponding to (84) is 

K = « + c x [0 - (2 - c x ) e] + c y [ 0 - (2 - c y ) e] 

+ cJO - (2 - c z ) e] > 0 (A.4.1) 

This is equivalent to 

( c x - l) 2 + ( c y - l) 2 + (c 2 - l) 2 > 2 (A. 4. 2) 

for ( c x , c y , c z < 1). The boundary of the positivity-preserving region is thus a portion of 
a sphere of radius sjl . For c x = c y = c z = c , this means that 

c < 1 - s/2/3 * 0.184 (A.4.3) 

compared with c < 1 - ^1/2 « 0.293 in the two-dimensional case. Clearly, (A. 4. 2) 
represents a rather severe restriction on the allowable time step. 


APPENDIX 5 

THREE-DIMENSIONAL LIMITER 

The three-dimensional generalisation of (130) and (131) can be derived in a straight-forward 
manner. For finding the slope 5), refer to Figure A. 1 . 1 . By analogy with (112), the update 
algorithm can be written 

<t>p = + C x C y + C y C z + ~ + C > + C z ) S l] 

+ C x( S l - C y - + C y ( S l ~ C Z ~ + C z( S i ~ C x ~ C y)<t> B 

+ C xCy <t>sw + c y c z <f> B + C Z C, <f> BW (A. 5.1) 

and this is required to be nonnegative. The $ P -term can be eliminated by choosing 

S = (‘ » ^ ^ C ,Q (A.5.2) 

‘ ( c , * c , * c ,) 

Note that this is consistent with the two-dimensional formula (c z = 0) and the one- 
dimensional case ( c y = c z = 0). The coefficient of <t> w , for example, then becomes 
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/I 2 2 v 

C A l ~ C y C z ~ C > ~ C z) 

(c x * c y + c z ) 

(A. 5. 3) 

with similar formulas for the coefficients of 0 S and 4> B obtained by permutation. These 
coefficients become negative in a rather significant portion of the ( c x , c y , c z ) space for which 
5i > 1, thus further restricting the PP region. 

The analogue of (126) becomes 

4>p = 0,(1 + c x c y + c y c z + c z c x - c x Sj - c y S l - c z Si) 


+ c x (S 2 X - c y - c z ) 4> w + c y (S 2 y - c z - c x ) <f> s + (Si 

- - c 7 ) 0 B 

+ C y $ SW + C Z 4* BS + C z C x $ BW 

(A.5.4) 

Eliminating 4> w , <f> s , and <f> B gives 


S 2 X = c y +c z 

(A.5.5) 

SI - + c. 

(A. 5. 6) 

Si = c x +c y 

(A.5.7) 

resulting in 


<t>; = 0,(1 - c,c y - c y c z - c z c x ) 


+ C x C y 0 SW + C z 0B5 + C z C x 

(A.5.8) 

representing a nonnegative weighting, provided 


c c + c c + c c <1 

x y y z z X 

(A.5.9) 

From the S 2 formulae it can be seen that, unfortunately, one of the S 2 s becomes unity 
whenever any two Courant number components are simultaneously equal to 0.5. For some 
applications this may be too restrictive. Further generalisations, re-introducing the non- 
symmetrical /function in order to find less restrictive limiter constraints, are currently being 


explored. 
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APPENDIX 6 

UPWINDING STRATEGY WITHOUT IF STATEMENTS 

Clearly, in order to take account of all combinations of signs of the Courant number 
components at any given face, the upwinding strategy could become quite complicated. 
Using IF statements to branch to a different formula for each combination is conceptually 
straight-forward but computationally very expensive. Figure A. 6. 1 shows the four combina- 
tions for a given face in the two-dimensional case, using a generic stencil notation 
representing upwind ([/), downwind (£>), central (C), etc., similar to Figure 17(a) for the 
one-dimensional case. Assume that the face in question is the left face of a particular control 
volume centred at node (I, J). The following upwinding strategy avoids IF statements, 
relying instead on the signs of the Courant number components at the face, and the 
FORTRAN nearest-integer function, NINT. 

Figure A. 6. 2 shows the labelling of node-values for the west face of a control volume 
(I, J) according to the signs of c xw and c yw . Note that any of the cases in Figure A. 6.1 can 
be viewed in the form of Figure A. 6. 2 (imagine a case in A. 6.1 drawn on a sheet of 
transparent film; by rotating and/or turning over the sheet, the configuration can always be 
made to conform to that in A. 6. 2). Using the nearest-integer function, NINT, the node 


labels can be described in terms of the signs of c xw and c as follows: 

IWU = I - NINT [0. 5 + SIGN(1.5, c x j] (A. 6.1) 

IWC = I - NINT[0.5 + SIGN(0.5, c w )] (A. 6.2) 

IWD = I - NINT [0. 5 - SIGN(0.5, c xw )] (A. 6.3) 

JWU = J - NINT[SIGN(1.0, c y j] (A. 6.4) 

and 

JWD = J + NINT[SIGN(1.0, c yw )] (A. 6.5) 

The six node-values of the generic third-order west-face stencil are given by 

<t>l D = <MIWC, JWD) (A.6.6) 

(t>u = <£(IWU, J) (A. 6.7) 


<H IWC, J) 


(A. 6.8) 
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<t>l = <^>(IWD, J) 


(A.6.9) 

<t>cu = *(IWC, JWU) 


(A.6.10) 

and 

<l> w D u = ^(IWD, JWU) 

The UTOPIA west face-value corresponding to (58) is then given by 


(A.6.11) 

Icl n — c 2 ) Icl 

4> = ^ - 1 GRADNW - V xw ' CURVNW - 1 yw ' 

LlN 2 6 2 

- GRADTW 

- |C ^' [(1 - 1 c yw 1 ) CURVTW - (1 - |c,J)TWISTW] 

(A.6.12) 

where 



* K) 


(A. 6. 13) 

GRADNL = 4>d ~ 4>c 


(A.6.14) 

CURVNL = </>£ - 2<t>c + <t>l 


(A. 6. 15) 

GRADTL = <t> w c - <i>cv 


(A. 6. 16) 

CURVTL = 4> w cd - 2<t> w c + 4> w cu 


(A.6.17) 

and 

TWISTL = 4>l ~ ~ 4>lu + <t>cu 


(A.6.18) 

Higher-order formulas corresponding to (59) can be designed in an obvious way. Note 
particularly the use of the absolute-values of the Courant number components in (A. 6. 12). 
A similar procedure is easily constructed for the other (south) face (in two dimensions). And 
it should be clear that the extension to three dimensions follows the same strategy for each 
of the three face-values involved for any given control volume. 
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Figure 1 . — Schematic of exact solution for 8 = 45°. The small cross (+) indicates the geometric centre of the domain at grid point (1 6, 1 6). Note 
that boundaries correspond to control -volume cell faces. Periodic boundary conditions imply: <J>(32, j) s 4>(1 , j), 4>(i, 32) s <j>(i, 1), etc. Positive 
"density" contours in this and all subsequent contour-plots are shown in increments of 0.1 . 
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(<|> = -0.001 , -0.05). 4> max = 0.928, <}> mjn = -0.091 . At = V2/1 24. Negative contours are: 4> = -0.0001 , -0.1 , -0.2, -0.3. 


Note the phase lag, slight overshoot (<J> max = 1 040), and large 
undershoots (<J> mjn = -0.366). 
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Figure 5. — Stability region (shaded) for the coordinate-wise 
first-order method. The scheme is stable for |c x | + |Cy| £ 1 . 



Figure 4. — Results for first-order upwinding used coordinate-wise for 
0 = 45°; c x = c y = 0.5; N = 62, At = >f2/62. 4> max = °- 473 *. no 
undershoots. 



Figure 6.— Advective flux at an angle into the west face of a 
control volume, suggesting dependence on the transverse 
gradient 
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Figure 7. — Results for first-order upwinding including GRADT terms; 
0 = 45°; c x = c y = 0.5; N = 62, At = V2/62. Note the good isotropy 
in comparison with Figure 4. 4> max = 0.366; no undershoot. 



Figure 8. — Calculation of the time-averaged west-face 
value assuming a linear (x, y) dependence in the vicinity 
of the face, collocated at the three nodes shown 
(for c x , Cy > 0). 
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Figure 9. — Results ol the generalised two-dimensional Lax-Wendroff 
scheme including GRADT terms; 0 = 45°; c x = c y = 0.25; N = 124, 

At = ''/2/1 24. = 0.867, <t> min = -0.087. 



Figure 10. — Six-point stencil used for deriving the two-dimen- 
sional Lax-Wendroff scheme by back-tracking along an 
advective characteristic. Note the upwind bias (for the case 
shown, c x , Cy > 0). 



Figure 1 1 . — Upwind-biased ten-point stencil used for the interpol- 
ation of 4> n (£, t\) given by Equation (44). 


-e-# 
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Figure 1 2. — Results for the uniformly third order polynomial interpola- 
tion algorithm (UTOPIA); O = 45°; c x = c y = 0.25; N = 1 24, At = V2/1 24. 
6 max = 0.885, 4> mln = -0.004. Note the good isotropy and phase 
accuracy. 



Figure 13.— UTOPIA results for 0 = tan-1(1/2) - 26.565°; c x = 0.5, 
Cy = 0.25; N = 1 24, At = VE/1 24. <t>max = 0.872, 4>mln = -0 006. 
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Figure 1 4. — UTOPIA results for 0 = 0°; c* = 0.5, c y = 0; N = 62, Figure 1 5. — Results using the one-dimensional QUICKEST scheme 

At = 1/62. <t>max = 0.958, 4> min = -0.002. coordinate-wise (Davis & Moore's method) for 8 = 45°; 

c x = c y = 0.25; N = 124, At = V2/124. <j>max = 1 - 167 . <l>min = -0.111. 







Figure 16. — ■Quasi-third-order" method consisting of QUICKEST 
together with GRADT terms (only) for 0 = 45°; c x = c y = 0.25; 

N = 1 24, At = V2/124. = 0.880, <)> m i n = -0.010. 




Figure 17. — Definition of geheric upwind (U), downwind (D), and central (C) node-values, depending on the 
sign of u. 
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(b) Exponential interpolation. 
Figure 18. — Interpolation through the same three (normalised) node-values. 


p wl 


T^ 


Figure 1 9. — Schematic one-dimensional control-volume show- 
ing east and west face-values. 



Figure 20. — One-dimensional universal limiter constraints 
corresponding to (63), (66), and (67). Note the slopes: 

Si = 1/c, and (for later reference) S 2 = 0. The case shown 
is for c = 0.25. 



Figure 21 . — Estimation of the time-averaged face- 
value, <f>j, assuming linear variation between 4> 
4>", and constant advecting velocity. 


03 
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Figure 22.— TVD flux-limiter constraints in terms of <j> " vs. Figure 23.— TVD flux-limiter constraints in terms of 4> f vs. 4>£. 

The case shown is for c = 0.25 (S«| = 1 .75, S 2 = 0.25). 



Figure 24.— Upper constraint slopes, Sj and S 2 , for ULTIMATE 
(heavy lines) and TVD schemes (light lines). 
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Figure 26.— Results using the one-dimensional TVD (CLAM) scheme 
coordinate-wise for 0 = 45°; c x = c y = 0.26; N = 1 24, At = V2/1 24. 
<t>max = 0.754, no undershoot 



Figure 27.— One-dimensional ULTIMATE QUICKEST scheme used 
coordinate-wise for 6 = 45°; c x =c y = 0.25; N = 1 24, At = V2/1 24. 


4* max = 0.985, = -0.1 24. 
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Figure 28 . — Normalised "front" passing obliquely through a control 
volume. 



Figure 29. — Restriction on the allowable component Courant numbers 
for positivity-preserving schemes using TVD limiters in each coor- 
dinate direction. Boundaries are quarter-circles with centres at 
(± 1 ,± 1 ). 



Figure 30. — Contours in the (c x , Cy) plane of the first 2D limiter slope, 
S-j, according to Equation 97. 
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Figure 31 . — Two-dimensional limiter of (97) applied to the QUICKEST 
scheme used coordinate-wise for 0 = 45°; c x = c y = 0.25; N = 1 24, 
At = V2/1 24. 4> max = 0.874, no undershoot 



Figure 32.— The two-dimensional limiter of (97) applied to UTOPIA 
for 0 = 45°; c x = c y = 0.25; N = 1 24, At = V2/1 24. 4> max = 0.777, no 
undershoot. 
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Figure 33. — Second two-dimensional limiter portrayed in the 
(<j>£ , iji®) plane, according to (119) and (120). The case 
shown is for c x = c y = 0.25 (S 1 = 2.125, S 2 = 0.25). 
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Figure 34. — Contour lines of S-j according to (130). Note The "flat" 
regions (S-j -> 1) towards the comers of the diagram. 
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Figure 35. — The two-dimensional limiter of (130), (131) applied to 
UTOPIA for 0 = 45°; c x = c y = 0.25; N = 1 24, At = V2/1 24. 

<|) rnax = 0.782, no undershoot. 



Figure 36. — The two-dimensional limiter of (130), (131) applied to 
UTOPIA for 0 = tan-1 (1/2) » 26.565°; c x = 0.25, c y = 0.125; N = 248, 
At = V5/248. 4> max = 0.752, no undershoot 
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Figure 37. — The two-dimensional limiter of (1 30), (1 31 ) applied to Figure 38. — The two-dimensional limiter of (1 30), (1 31 ) applied to 

UTOPIAfor 0 = tan-1 (1/3) * 18.435°; c x = 0.25 = 3 c^,; N = 372, UTOPIA for 9 = 0°; c x = 0.5, c y = 0; N = 62, At = 1/62. 

At = ViO/372. = 0.731 , no undershoot 4>max = 0.927, no undershoot 
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Figure 39.— The two-dimensional limiter of (130), (131) applied to 
UTOPIA for 6 = 45°; c x = c y = 0.5; N = 62, At = V2/62. 

<J> max = 0.715, no undershoot. 
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adjacent faces. 
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(b) Three outflows. 
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(c) Three inflows. (d) Inflows (outflows) on 

opposite faces. 

Figure 40. — Four inflow-outflow configurations. 





Figure A.1 .1 . — Twenty-node three-dimensional molecule for third 
order upwinding. For the case shown, c Xl Cy, c z >0. 


Figure A.1 .2. — Upwind-biased ten-node molecule used for 
computing the west-face advected value. The case shown 
is for c x , Cy, c z >0. 



Figure A2.1. — Variation of the normal gradient at the west face due 
to curvature in the advected profile. 
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Log-jo (CPU seconds) 


Figure A.3.1. — Computational efficiency diagram (error versus cost, on a log-log scale) for the two-dimensional generalised 
Lax-Wendroff scheme, the one-dimensional QUICKEST scheme used coordinate-wise together with transverse-gradient terms, 
and the uniformly third-order polynomial interpolation algorithm (UTOPIA). 




(a) c x > 0, c y > 0. (b) c x < 0, c y > 0. 



(C) C X > 0, Cy < 0. 



(d) c x < 0, c y < 0. 


4>cd 
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Figure A.6.1 . — Four possible two-dimensional stencils depending on Figure A.6.2. — Generic two-dimensional stencil for the west 

the signs of the velocity components at the face in question. face of the control volume. 
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